This material is based on a document created by Claudia Engle, available at https://cengel.github.io/rspatial/4_Mapping.nb.html .
Lovelace, Nowosad, and Muenchow have a nice open source book on Geocomputation in R. https://geocompr.robinlovelace.net/
It covers much more that just visualization, but Chapter 8 covers a bit about making maps. https://geocompr.robinlovelace.net/adv-map.html
Download the HomeworkMapsPt2.nb.html file from
Moodle
From the notebook file, , create the associated .Rmd file and open in RStudio.
Replace the “Your Name Here” text in the author:
field with your own name.
Supply your solutions to the homework by editing
HomeworkMapsPt2.Rmd.
When you have completed the homework and have
checked that your code both runs in the Console and
knits correctly when you click Knit HTML, rename the R
Markdown file to HomeworkMapsPt2_YourNameHere.Rmd, and
submit both the .Rmd file and the .html output
file on Blackboard. (YourNameHere should be changed to your own
name.)
| Keystroke | Description |
|---|---|
<tab> |
Autocompletes commands and filenames, and lists arguments for functions. |
<up> |
Cycles through previous commands in the console prompt |
<ctrl-up> |
Lists history of previous commands matching an unfinished one |
<ctrl-enter> |
Runs current line from source window to Console. Good for trying things out ideas from a source file. |
<ESC> |
Aborts an unfinished command and get out of the + prompt |
Note: Shown above are the Windows/Linux keys. For
Mac OS X, the <ctrl> key should be substituted with
the <command> (⌘) key.
Instead of sending code line-by-line with
<ctrl-enter>, you can send entire code chunks, and
even run all of the code chunks in your .Rmd file. Look under the
Run your code in the Console and Knit HTML frequently to check for errors.
You may find it easier to solve a problem by interacting only
with the Console at first, or by creating a separate .R
source file that contains only R code and no Markdown.
As this is a new day, and a new file, we still have some preliminary setup that is required.
So … yesterday ggplot2 version 3.0 becamee
available, which includes the geom_sf support. It also
incorporates viridis as its default colormap for ordered
variables.
You should update that package from inside of Rstudio.
I tried to see if ggmap had incorporated a fix, and it had not, but … I won’t be using that stuff today. On the other hand, I generally load ggplot2 out of habit when I am working in R.
We will give preference to the “tidyverse” of packages related to “Tidy Data.”
Pt 1 loaded several spatial packages. For this part, I am commenting out several of those package loads so that you better understand the particular dependencies.
## additional libraries needed for R code examples
library(sp)
library(sf)Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE
#library(rgdal)
#library(classInt)
library(RColorBrewer) # not a strictly a spatial ligrary
library(viridis) #my current "go to" color scalesLoading required package: viridisLite
#library(ggmap)
library(leaflet)
library(tmap)
#library(GISTools)Some general tools that I almost always load:
library(knitr) # to create nice documents in R
library(tidyverse) # loads ggplot2, dplyr,tidyr,readr,purr,tibbleRegistered S3 methods overwritten by 'dbplyr':
method from
print.tbl_lazy
print.tbl_sql
── Attaching packages ────────────────────────────────────────────────────────── tidyverse 1.3.2 ──✔ ggplot2 3.3.6 ✔ purrr 0.3.4
✔ tibble 3.1.8 ✔ dplyr 1.0.10
✔ tidyr 1.2.1 ✔ stringr 1.4.1
✔ readr 2.1.2 ✔ forcats 0.5.2 ── Conflicts ───────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
library(broom) # because I find it useful
library(forcats) # working with categorical variablesWe will continue to play using the Philly data, so that we maintain a sense of continuity.
Philly3.Make sure that you place an unzipped copy of that directory as a subdirectory in your working directory.
For this lecture we, we will focus ourselves entirely on the newer,
sf objects, so I won’t even load the older format.
Recall the meaning of the “added” variables:
philly_sf <- st_read("Philly3/Philly3.shp")Reading layer `Philly3' from data source
`C:\Users\rashi\Downloads\CU\Fall\info.Visualization\map_hw\Mapping files-20221129\Philly3\Philly3.shp'
using driver `ESRI Shapefile'
Simple feature collection with 381 features and 12 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 1739498 ymin: 457288.9 xmax: 1764018 ymax: 490539.6
Projected CRS: Albers
plot(philly_sf)Warning: plotting the first 10 out of 12 attributes; use max.plot = 12 to plot all
sfis an important package for data manipulation and integration with other programs (hence - conformance to formal ISO standards). If you are thinking that you need to manipulate the underlying data (or ask geographic questions), it is probably the tool. The other packages we are using focus on visualization.
philly_sf?tmapThe tmap was specifically designed for making thematic
maps, where the focus is on geospatial map as a context to understanding
some other data (such as we might be doing with this homicide data) It
borrows from the ggplot syntax (grammar of graphics) and is specifically
designed to make creation of thematic maps more convenient. It takes
care of a lot of the styling and aesthetics. This reduces our amount of
code significantly. We only need:
tm_shape() and provide the sf object as
argument directly. This is followed by+, andtm_polygons where we set
tmap does not require an
sfobject, but it can certainly use them.
Check the help tm_shape() and tm_polygons
to see how you can set optional parameters (map, break, title) and try
on your own before you look at my code.
tm_shape(philly_sf) +
tm_polygons("HOMIC_R", style="quantile", title="Philadelphia \nhomicide rate \nper 100,000")The “standard” for posting packages to CRAN is that they should include a vignette - which describes a bit about using the package. Take a look at https://cran.r-project.org/web/packages/tmap/vignettes/tmap-nutshell.html . It can give you an idea of some additional capabilities.
tm_shape(philly_sf) +
tm_polygons("mdHHnc_", style="quantile", title="Median Household Income")spplot(philly,"mdHHnc_")tm_shape(philly_sf) +
tm_polygons(c("HOMIC_R","mdHHnc_"), style="quantile")Notice the challenge of trying to understand the relationship between these two variables, even when plotted side by side.
In the vignette, it points out that there are other
marks that one can use (besides the polygons). Here is a
bubble representation of the income:
tm_shape(philly_sf) +
tm_bubbles("mdHHnc_",perceptual=TRUE)tm_shape(philly_sf) +
tm_polygons("HOMIC_R",style="quantile")+tm_bubbles("mdHHnc_",perceptual=TRUE)tm_shape(philly_sf) +
tm_bubbles("mdHHnc_",perceptual=TRUE)+
tm_polygons("HOMIC_R", style="quantile",alpha=.4)There is a second vignette for tmap (https://cran.r-project.org/web/packages/tmap/vignettes/tmap-modes.html) which describes an alternative drawing mode. Above, we have been using a representation suitable for static maps. But tmap can also access leaflet tools to draw html widget based maps that allow for a bit more iteractivity.
We can easily switch from “plot” mode into “view” mode and call the last plot, like so:
tmap_mode("view")
tmap_last()Cool huh?
The tmap library also includes functions for simple
spatial operations, geocoding and reverse geocoding using OSM. For more
check vignette("tmap-nutshell").
a=tmaptools::geocode_OSM("Rocky statue, Philadelphia, PA",as.sf=T)
tm_shape(philly_sf) +
tm_polygons("HOMIC_R", style="quantile")+
tm_shape(a) + tm_text("query",size=1)leafletLastly, leaflet[^7] makes use of the widely known ‘Leaflet’ JavaScript library, “the
leading open-source JavaScript library for mobile-friendly interactive
maps”. We have already seen a simple use of leaflet in the
tmap example.
The good news is that the leaflet library gives us loads
of options to customize the web look and feel of the map.
The bad news is that the leaflet library gives us loads
of options to customize the web look and feel of the map.
You don’t have to, but it makes the code more accessible if you use
the pipe operator
%>% to chain the elements together when building up
a map with leaflet.
Let’s build up the map step by step.
We have spoken of coordinate reference frames before. Leaflet requires the spatial data to have preassign reference and likes the WGS84 data. So as first step, I will transform our cordinates into that crs.
ph_WGS = philly_sf %>% st_transform(4326)We can visualize,
# first try... ops what happened here
leaflet(ph_WGS) %>%
addPolygons()but let’s tune what we want for the polygons:
Map the homicide rate. For this we provide several parameters to the
addPolygons() function that:
- remove stroke (polygon borders)
- set a fillColor for each polygon based on `HOMIC_R` and make it look nice by adjusting fillOpacity and smoothFactor (how much to simplify the polyline on each zoom level). The fill color is generated using the `colorQuantile()` function, which takes the color scheme and the desired number of classes. `colorQuantile()` then returns a function that we supply to `addPolygons()` with the name of the value we want to map to constuct the color scheme.
- add a popup with the `HOMIC_R` values. We will create as a vector of strings, that we then supply to `addPolygons()`.
pal_fun <- colorQuantile("YlOrRd", NULL, n = 5) # helper function to allow coloring by variable value
p_popup <- paste0("<strong>Homicide Rate: </strong>", ph_WGS$HOMIC_R)
leaflet(ph_WGS) %>%
addPolygons(
stroke = FALSE, # remove polygon borders
fillColor = ~pal_fun(HOMIC_R), # set fill color with function from above and value
fillOpacity = 0.8, smoothFactor = 0.5, # make it nicer
popup = p_popup) # add popupAdd the default basemap, which is OSM, with addTiles().
This you can do!
leaflet(ph_WGS) %>%
addPolygons(
stroke = FALSE,
fillColor = ~pal_fun(HOMIC_R),
fillOpacity = 0.8, smoothFactor = 0.5,
popup = p_popup) %>%
addTiles()Finally, I have added for you a control to switch to another basemap and turn the philly polygon off and on. Take a look at the changes in the code below.
leaflet(ph_WGS) %>%
addPolygons(
stroke = FALSE,
fillColor = ~pal_fun(HOMIC_R),
fillOpacity = 0.8, smoothFactor = 0.5,
popup = p_popup,
group = "philly") %>%
addTiles(group = "OSM") %>%
addProviderTiles("CartoDB.DarkMatter", group = "Carto") %>%
addProviderTiles("Esri.WorldImagery", group = "Sat") %>%
addLayersControl(baseGroups = c("OSM", "Carto","Sat"),
overlayGroups = c("philly")) So Leaflet ultimately has more capability, but a bit more work to get there.